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Majorana fermions are predicted to play a crucial role in condensed matter realizations of topo- 
logical quantum computation. These heretofore undiscovered quasiparticles have been predicted to 
exist at the cores of vortex excitations in topological superconductors and in heterostructures of 
superconductors and materials with strong spin-orbit coupling. In this work we examine topological 
insulators with bulk s-wave superconductivity in the presence of a vortex-lattice generated by a 
perpendicular magnetic field. Using self-consistent Bogoliubov-de Gennes, calculations we confirm 
that beyond the semi-classical, weak-pairing limit that the Majorana vortex states appear as the 
chemical potential is tuned from either side of the band edge so long as the density of states is 
sufficient for superconductivity to form. Further, we demonstrate that the previously predicted vor- 
tex phase transition survives beyond the semi-classical limit. At chemical potential values smaller 
than the critical chemical potential, the vortex lattice modes hybridize within the top and bottom 
surfaces giving rise to a dispersive low-energy mid-gap band. As the chemical potential is increased, 
the Majorana states become more localized within a single surface but spread into the bulk to- 
ward the opposite surface. Eventually, when the chemical potential is sufficiently high in the bulk 
bands, the Majorana modes can tunnel between surfaces and eventually a critical point is reached at 
which modes on opposite surfaces can freely tunnel and annihilate leading to the topological phase 
transition previously studied in the work of Hosur et al} . 

PACS numbers: 71.10.Fd, 75.10.Jm, 71.10.Pm, 75.40. Mg 



INTRODUCTION 



Majorana fermions, quasi-particle excitations which 
are their own antiparticle, were originally proposed in 
high-energy physics but^ have now arrived at the fore- 
front of condensed matter physics where they serve as 
non-Abelian anyons which form the backbone of topolog- 
ical quantum computing architectures^"—. Within con- 
densed matter physics, there exist many candidate sys- 
tems which are predicted to harbor Majorana fermions. 
One of the earliest of such candidates is the fractional 
quantum Hall effect at filling factor v = the physics 
of which may be described by the Moore-Read pfaf- 
fian wavefunction^^. While this state is yet to be ex- 
perimentally confirmed, tantalizing evidence observed in 
tunnehng in quantum constrictions points to the fact 
that the ^ = f fractional quantum Hall state does 
possess non-Abelian statistics as would be necessitated 
by the presence of Majorana fermions^^. Beyond the 
fractional quantum Hall states, other possible systems 
thought to contain Majorana fermions are the px + ipy 
superconductors^^i^ T 1 1 1 13^ 14 ^j^g^.^ relevant Majorana 
modes are predicted to appear as bound-states on ex- 
otic half-quantum vortices, which were recently observed 
in magnetic force microscopy experiments performed in 
Sr2Ru04^^. In addition to fractional quantum Hall states 
and superconductors, there has been an abundance of 
proposals to realize Majorana fermions in materials with 
strong-spin orbit coupling. Notable examples are prox- 
imity induced superconductivity in 3D topological in- 
sulators (TIs)^^, bulk superconductivity in doped TI^ 



and semiconductors coupled proximity coupled to s-wave 
superconductors^—. Indeed, the latter proposals have 
led to exciting measurements in high mobility quantum 
wire - s-wave superconductor systems^. 

In this article we will focus on the two mechanisms 
proposed in TI materials. As is now well-known, TIs 
are materials which possess an insulating bulk but con- 
tain robust metallic states that are localized on their 
surfaces^^^^. We will consider time-reversal invariant 
3D topological insulators which harbor an odd number 
of massless Dirac cones on each surface. As mentioned, 
there currently exist two proposals that utilize topologi- 
cal insulators as a platform for the observation of Majo- 
rana fermions. The first of which is to consider an 5- wave 
superconductor /topological insulator heterostructure in 
which a superconductor is coupled to the topological in- 
sulator via the proximity effect and subjected to a vortex- 
producing magnetic field. In Fu and Kane's pioneering 
work,^^ they show that in the 5- wave superconductor /TI 
heterostructure, the interface between the TI and the su- 
perconductor behaves similar to a spinless chiral p-wave 
superconductor yet without breaking time reversal sym- 
metry. As such, Majorana fermions will reside in the 
vortex core^i^ so long as the quantized magnetic flux 
lines penetrating the system are broad^^. 

On the other hand, another strategy to realize Majo- 
rana fermions is to consider vortex bound-states in a 3D 
TI with bulk 5-wave superconductivity^^i^J^. Bulk su- 
perconductivity in doped topological insulators has been 
observed in recent experiments that dope Bi2Se3 with 
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the order parameter is stih under debate, but the two 
most probable options are s-wave, or an inter-orbital 
topological pairing parameter^^. There is an opportu- 
nity to observe Majorana fermions in both cases, but for 
the purpose of this work we will only consider the s-wave 
case. Recent work^ reveals that, while doped topological 
insulators that develop 5- wave pairing may harbor Majo- 
rana bound states in the vortices, the Majorana fermions 
do not survive for all doping levels. Specifically, there 
exists a critical chemical potential, /ic? at which point 
the system undergoes a topological (vortex) phase transi- 
tion. This phase transition can be regarded as a topology 
change in the ID electronic structure of vortex lines from 
a system which supports gapless end states to one that 
does not. Therefore, only at chemical potentials below 
the critical value, the doped superconducting Bi2Se3 sup- 
ports Majorana modes at the vortex ends (places where 
vortex lines intersect the surface). The work of Ref. Jj 
provided a semiclassical treatment in the infinitesimal 
pairing limit, but such an approach might be inadequate 
to capture important quantum effects. One such effect 
that cannot be determined this way is the zero-point en- 
ergy contribution of the vortex core states which could 
shift the states away from the zero energy gapless point 
and invalidate the previous analysis. In the weak-pairing 
limit the physics is determined by the structure exactly at 
the Fermi-surface, and it is possible that the energy spec- 
trum away from the Fermi-level can serve to renormalize 
the location of the critical point. If the critical point is 
sufficiently shifted toward the band-edge, it could be that 
there is never a viable doping-range over which the Ma- 
jorana fermions can be observed. Calculations provided 
in this article go beyond the semi-classical, and infinites- 
imal weak-pairing limit are thus essential to confirm the 
previous results. 

In both TI based approaches we have thus far dis- 
cussed, there is an additional assumption underlying the 
resultant physical predictions, namely that the vortices 
are completely isolated. However, this may not be the 
most appropriate or experimentally relevant picture for 
the observation of Majorana states in either type of TI 
approach. Thus, in this article, we examine the behavior 
of 3D topological insulators with bulk 5-wave supercon- 
ductivity in the vortex lattice limit as a function of dop- 
ing level. The paper is organized in the following fashion: 
In Sections [Til and HlH we introduce the 3D topological 
insulator model Hamiltonian which is used for each of the 
subsequent calculations, and the background for the self- 
consistent calculations respectively. In Section IIV( we 
present the results of our calculations for three separate 
geometries (a)periodic boundary conditions with vortex 
rings (b) open boundary conditions with vortex lines 
terminating on the TI surface (c) an inhomogeneously 
doped heterostructure with open boundary conditions. 
We find that as the chemical potential is adiabatically 
moved from the gap into the bulk bands, the Majorana 
states form when the density of states reaches is large 
enough to support a well- formed superconducting gap. 



As the chemical potential moves past this onset value we 
find that the vortices are localized on the surfaces but 
hybridize with neighboring vortices on the same surface 
giving rise to a dispersive low-energy quasi-particle spec- 
trum. As the chemical is pushed further into the bands, 
we find a critical chemical, /i^, at which inter-surface tun- 
neling is enabled through a gapless channel on the vortex 
line. The value is renormalized from that stated in Ref. 
1, but we find that even for strong attractive interac- 
tions that fit remains at a finite value of the bulk-doping. 
After the chemical potential exceeds /it, we find that a 
gap opens in the spectrum and there are no longer any 
low-energy localized modes remaining. Additionally, in 
Section [Vl we evaluate the superconducting gap equation 
in order to determine the relevant temperature scale on 
which these effects may be observed. Finally, in Section 
I VII we summarize our findings and conclude. 

II. MODEL HAMILTONIAN OF A 3D 
TOPOLOGICAL INSULATOR 

We will use a minimal, four-band Dirac-type model 
which, with the proper choice of parameter values, cap- 
tures the bulk, low-energy physics of known TI materials 
such as Bi2Se3i^i^: 

if„ = Mr», h, = J:!^1^^. (2) 
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where = (cA,t,r CA,i,r CB,t,r CB,i,rY '^^ ^ Sub- 
component spinor with A/B and t / i labeling orbital 
and physical spin respectively so that cj^ ^ ^ is the cre- 
ation operator for an electron with spin cr in orbital a 
at position r , J = ±af , ±a^, ±ai are vectors that con- 
nect nearest neighbors on a simple cubic lattice with lat- 
tice constant a, the vector F = T^x + F^^ + T^z with 
F" = (g) cr" and F^ = I, where a = x, 
and cr" are 2x2 Pauli matrices acting on orbital and 
spin degrees of freedom, respectively. We also define I 
as the 2x2 identity matrix and M = m — 36/a^ as the 
mass parameter which controls the magnitude of the bulk 
band gap. The Tl/trivial insulator phase depends on the 
chosen values for the parameters m and h and the TI 
phase has m/b > while the trivial phase has m/b < 0. 
By tuning the material parameters 7, 6, m, and a in Eq. 
([1]) one can model the low energy effective model for the 
common binary TI materials^^^^^. Since we only address 
the qualitative effects stemming from the TI phase we 
will fix the parameters to be 6 = a^(leV), 7 = a(leV) 
and m = 1.5 eV in terms of the lattice constant a so that 
M = —1.5 eV thereby ensuring that we are in the TI 
phase. 

With translation invariance and periodic boundary 
conditions in all x, y and z directions, it is often more 
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convenient to work in momentum space. In this case, 
we expect to see no gapless states due to the lack of a 
boundary. The Fourier transformed Dirac Hamiltonian 



in momentum space may be written as 

k 

where Ho{k) in Eq. (|3]) is a 4 x 4 matrix which we may 
write as 



Ho{k) = 





sin kz 







sin kz 



\ sin kx -\- i sin ky 



M + ^(^) sin kx -\- i sin ky 
sin kx — i sin ky —[M-\-g{k)] 
— sin kz 



sin kx—i sin ky \ 
— sin kz 


-[M + ^(^)] J 



(4) 



with g{k) = cos kx + cos ky + cos /c^. If we instead choose 
open boundaries along one direction there will be robust 
gapless edge states on those boundary surfaces for the 
same choice of model parameters. 



where = , ^t)"^ is now an 8-component Nambu 
spinor, and A(r) denotes a 4 x 4 pairing matrix. In this 
expression, the interaction —\U\ has been absorbed into 
the pairing matrix A(r), which we write as 



III. TOPOLOGICAL INSULATORS WITH 
BULK 5^- WAVE SUPERCONDUCTIVITY 

In this work, we are interested in the properties of 
doped TIs which become intrinsically superconducting at 
low temperature. In a doped topological insulator, like 
any other metal, when the chemical potential is in the 
conduction or valence band an attractive interaction will 
lead to the formation of superconductivity and generate 
a superconducting gap at the Fermi surface. In order to 
study the formation of superconductivity in doped TI we 
add an attractive Hubbard-type density-density interac- 
tion to the Hamiltonian in Eq. ([T]) : 



Hint = -\U\^n^,^i^f^ 



(5) 



^A,a,f^A,a,f^^B,a,f^B,a,f parameter 
— \U\ represents the attractive intra-orbital interaction. 

At the mean-field level the interaction term may be 
decoupled as^: 
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where Ac^^^ = {ca,i,rCa,^,r) is the standard intra-orbital 
s-wave pairing order parameter. Combining this with 
Eq. ([1]), we get the Bogoliubov-de Gennes (BdG) Hamil- 
tonian: 
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To study the bulk superconductivity we will assume /j.^ = 
/i is uniform throughout the material for simplicity. 

The BdG Hamiltonian of Eq. ([6]) can be diagonalized 
by applying a Bogoliubov transformation as — 



E 



Vr),_r 



(8) 



where n labels the eigenstate index. Plugging the trans- 
formation into Eq. (|6]), we have 



HBdG/Z 

n 

= E 



Vri.r 



En 







-En 



Vn.r 



(9) 



This indicates that the eigenvectors associated with En 
{—En) of the above BdG equations are {un,r^Vn,r)^ 
[(— 'U* ii* _.)^]. The mean-field pairing order parame- 
ters are obtained via 



= Un^r^n.f^^^^ 



PEn 



(10) 



where j3 = l/ksT. Once the pairing order parameter 
are determined initially it is plugged back into the BdG 
Hamiltonian given in Eq. (|6|) and then HBdG is diago- 
nalized again as shown in Eq. (j9j). The process continues 
until we reach self-consistency and we have a convergent 
Aq, for all r. We note that in our numerical calculations 
we use a small, non-zero temperature in order to avoid 
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FIG. 1. (color online). Intra-orbital pairing order parame- 
ters Aa (solid symbols) and Ab (hollow symbols) vs at 
different \U\. and \U\ are in units of eV. The mass term 
M = —1.5 eV. The system contains periodic boundary condi- 
tions in X, y and z directions. The simulations are performed 
on a lattice grid of size 80a x 80a x 10a. 



divergences but this temperature is much smaller than 
the superconducting gap so as not to affect the physical 
results. 

In Fig. [1] we show the self-consistently determined 
intra-orbital pairing order parameter in the bulk as a 
function of at different |[/|, where, due to translation 
invariance, A^^^^ = Ac^. In this paper we will only con- 
sider p-doping (/i < 0), but the particle- hole symmetry 
of the model Hamiltonian in Eq. ([1]) makes the electron- 
doped case similar in nature. With M = —1.5 eV in 
Eq. ([T]), the top of the bulk valence band is located at 
/jjy = —0.5 eV and the total size of the insulating gap is 
1.0 eV. We see from Fig. [T]that when the chemical po- 
tential is in the gap where there is no carrier density with 
which to form Cooper pairs and the resulting pairing po- 
tential is zero. When the chemical potential enters the 
valence band a Fermi-surface develops, and low-energy 
states become available to pair. However, when the den- 
sity of states at the Fermi- level is insufficient, the size 
of the pairing potential will continue to be exponentially 
small. As we see in Fig. (TJ a significant pairing potential 
does not form until is well above the valence band 
edge, 

This result matches standard BCS phenomenology and 
represents the point of inception for the remainder of the 
paper. To be specific, in Ref. T Hosur et al. used a semi- 
classical treatment to show that a vortex in the supercon- 
ducting phase of a doped topological insulator exhibits a 
topological phase transition as the chemical potential is 
tuned through a critical value. The two phases separated 
by this transition are gapped and differ by the presence or 
absence of Majorana modes at the ends of the vortex i.e. 
where the vortex line intersects the surface of the TI. At 
the transition point the vortex line becomes gapless and 
provides a channel which allows the Majorana modes to 
annihilate one another by tunneling in-between the op- 
posing surfaces. In their treatment, however, there is an 



assumption of adiabaticity as it is always assumed that 
at any chemical potential other than the critical chemi- 
cal potential, there is no gapless channel to hybridize the 
Majorana mode. This seems innocuous, but one has to 
remember that the arguments rely on the adiabatic con- 
nection between a gapped insulating phase and a gapped 
superconducting phase. The assumption enters when 
one considers the behavior of the system as the chem- 
ical leaves the insulating gap and enters the bulk bands. 
Although this is a reasonable assumption within which to 
theoretically study the vortex phase transition, one may 
then ask what happens in the region where the chemical 
potential is not large enough to form a significant pairing 
potential, and there is finite density of gapless modes in 
the bulk. In other words, how does the Majorana mode 
emerge out of the bulk gapless states? This question is 
certainly relevant for experiments where a finite size TI 
sample is used. Our self-consistent solution of the BdG 
equations in the vortex lattice can present a clearer pic- 
ture of the appearance of the Majorana modes and the 
vortex phase transition than the previous semi-classical 
analysis. 



IV. VORTEX LATTICES IN 
SUPERCONDUCTING PHASE OF DOPED 
TOPOLOGICAL INSULATORS 

The self-consistent BdG formalism is in a real space ba- 
sis and thus can be also used to study the vortices in the 
superconducting phase where the order parameter will be 
non-uniform. To induce vortices, we consider the system 
under a uniform magnetic field B = Bz. When electrons 
are hopping on the x?/-plane this generates a Peierls phase 
factor, and the BdG Hamiltonian becomes^^'"^ 



rr f A(r) 



t / Hse-'^^^ 



r,6 



(11) 



where r]^ denotes the extra phase given by the vector 
potential A{r) induced by the magnetic field through B = 
V X A{r): 



Vr- 



TT 
$0 



r+5 



Af^^ • dr 



(12) 



where is the superconducting flux quantum; $o = ^• 
In the following discussion, we choose the Landau gauge, 
i.e. A{r) = {A^,Ay) = {0,Bx). 

We will treat the system as a type-II superconductor 
in a vortex lattice state. In each magnetic unit cell, the 
amount of magnetic flux is 2^0, so that each unit cell 
carries two superconducting vortices^. We designate the 
size of each magnetic unit cell as Ixa x lya x l^a using 
the integers k to denote the number of lattice sites in 
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FIG. 2. A 4 X 4 vortex lattice. In this example, the number 
of magnetic unit cells is Nx x Ny =4x2. Each black solid 
circle denotes a vortex location and each magnetic unit cell 
contains two vortices. 



each spatial direction. For our choice of geometry we 
will use square vortex lattices, and fix = lyf^- The 
corresponding magnetic field magnitude is 



B 



2^0 



(13) 



From this relation we can observe that stronger mag- 
netic fields bring smaller magnetic unit cells, in which 
vortices are closer each other. Therefore, the dilute vor- 
tex limit comes from applying very weak magnetic fields. 
As is standard for lattice calculations with uniform field, 
in order to see experimentally reasonable field sizes we 
would need to use a very large number of unit cells as, 
for example, the case when = ly = ^ gives a magnetic 
field on the order of thousands of Tesla. For our system 
sizes we have an un-physically large magnetic field on the 
order of 10^ T assuming a lattice constant of lA. This, 
however, will not affect the qualitative physics in which 
we are interested and we will not worry about this issue 
any further. 

We choose the entire system size as LxCl x Lya x l^a 
such that there are Nx x Ny magnetic unit cells, where 
Nx = Lx/lx and Ny = Ly/ly and A^^^ = 2A^^. Since each 
magnetic unit cell carries two vortices, the Lxa x Lya x /^a 
vortex lattice contains 2NxNy vortices. In Fig. [2l we 
show a schematic illustration of a 4 x 4 square vortex lat- 
tice. By tuning sizes of the magnetic unit cells, we can 
study the vortex lattice at different external magnetic 
fields. In this paper, we set Ix > S to avoid strong over- 
lap between vortices but Ix < 12 due to computational 
limitations. 

We consider a system with periodic boundary condi- 
tions along the x and y directions. Although the vor- 
tices break lattice translation invar iance, we still have 
magnetic periodic boundary conditions for vortex lat- 
tices. In addition to the phases given by vector potentials 
e*^^, the magnetic periodic boundary conditions also con- 
tribute another phase factor when electrons are hopping 
across unit cell boundaries^^^^^. Suppose that in a 2D 
vortex lattice, the translation vector in units of a is writ- 
ten as R = {Xlx^Yly)^ where X = 0, • • • , A/'^^ — 1 and 
Y = 0^ - ■ ■ , Ny — 1 are integers. The coordinate of an 
arbitrary lattice site can be expressed as r + R, where 



r = (x, y) denotes the coordinate in units of the lattice 
site, a, within a magnetic unit cell, i.e. 1 < x < Ix and 
^ ^ y ^ ly Under the magnetic periodic boundary con- 
ditions, we can define the relation of the magnetic Bloch 
wave functions^ which have a periodic structure written 



Vn{r^ IxX) 
Un{r^ lyy) 

Vn{r^ lyy) 



27ri 



'vn{r) 



Un{r) 
Vn{r) 



(14) 



Here kx = and ky 



represent the x and y com- 
ponents of the magnetic Bloch wavevector. The phases 
6*^"= and e^^y arise from hopping to neighboring cells. 
Additionally, e is provided by the magnetic pe- 

riodic boundary conditions, or quasi-periodic boundary 
conditions^''. The BdG eigenstates (un^r^Vn^r)^ satisfy 
magnetic translation invariance under Eq. (p!4|) . 

The on-site pairing potential can be expressed as 
Aq,^^ = \Aa^f\e'^^^^'^ with a phase e*"^^^^ and amplitude 
|Ac^^f?|. In the presence of vortices, both the pairing po- 
tential and the phase are site-dependent. The supercon- 
ducting order parameters are suppressed near the vortex 
cores, and are restored to the bulk values away from the 
vortex cores. The spatial form of the pairing order pa- 
rameters Aa^r are determined self-consistently. We con- 
sider different attractive-Hubbard interaction strengths 
\U\ and uniform doping- levels distributed through the 
entire bulk. 

We study two different geometries for the vortex lat- 
tice. First we consider vortices oriented in the z-direction 
(along the applied magnetic field) with periodic bound- 
ary conditions along the x^y^z directions which yields 
vortex rings looping around the z-direction. In this ge- 
ometry we study the vortex phase transition where the 
vortex modes become gapless along the vortex rings. The 
second geometry we consider has open boundaries in the 
z-direction. In this case the vortex lines terminate at 
the open surfaces perpendicular to the z-axis and we can 
study the Major ana modes that can appear at the vortex 
ends in the topological phase. We compare these results, 
which are neither in the semi-classical, or infinitesimally 
weak-pairing limits, to the results studied in Ref. Jj which 
are in these limits. 



A. Periodic Vortex Rings in Vortex Lattices 

With periodic boundary conditions in all spatial direc- 
tions, we cannot directly study the Majorana modes that 
might appear at the vortex ends. However, we can indi- 
rectly study them by identifying the vortex phase transi- 
tion through a study of the low-energy modes along the 
vortex lines. As the chemical potential is tuned deeper 
into the band, the point where one of these modes be- 
comes gapless signals the location of a critical point. For 
this geometry the magnetic unit cell sizes we use are 
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FIG. 3. (color online). The energy spectra of the low-energy 
states vs for periodic boundary conditions along the z-axis 
(vortex rings) for (a) = 2 eV, (b) \U\ = 2.5 eV, and (c) 
\IJ\ = 2.8 eV. The systems have a A^^ x AT^ = 10 x 5 (100 
vortices) vortex lattice and the size of the magnetic unit cell 
is X L X Iz ^ 12a x 24a x 10a. 



IxXlyXlz = 12 x 24 x 10. We choose NxXNy = 10x5 unit 
cells so that there are 100 vortices in the vortex lattice. 

In Fig. [3l we present the evolutions of the low energy 
states vs for different interaction strengths \U\. We 
can identify two distinctly different doping regimes. In 
the first regime, the chemical potential lies in the valence 
band but below a value we call |/io| which signals the on- 
set of a well-formed superconducting gap discerned from 
our numerics. It should be noted that |/io| has no real 
intrinsic meaning (as it is strongly finite-size dependent) 
and only serves to indicate a common feature shared by 
all of our spectrum plots. Clearly, before the chemical 
potential hits the top of valence band there is no density 



of states to generate the superconducting gap, and all 
the states are gapped by the bulk insulating band gap. 
After the chemical potential hits the top of valence band, 
the superconducting pairing starts to form but it is ex- 
ponentially small in magnitude. Comparing the pairing 
strengths without vortices {i.e. Fig. [1]) to the vortex 
lattice case in Fig. [3l our numerics show that in the 
vortex lattice the pairing is more poorly formed over a 
larger range of doping. That is, the states at the Fermi 
level remain gapless with no superconducting gap forma- 
tion. Note that /io decreases with increasing \U\ which 
indicates that this point is sensitive to the point where 
the exponentially suppressed superconducting gap would 
turn on. In this regime, any localized Majorana modes or 
low-energy vortex core states are difficult to distinguish 
from the extended gapless metahic states in the bulk. 
The details of this regime are dominated by strong finite- 
size effects. One hinderance is that for cases where only 
a tiny pairing potential would form it is numerically chal- 
lenging to generate a convergent, self-consistent solution 
with vortices present. In the thermodynamic limit, we 
would expect to see a non-zero but exponentially small 
pairing gap as soon as the chemical potential hits the 
valence band. Here the situation is not so clear, and 
unfortunately, due to the computational limitations, we 
cannot glean a great deal of physical information from 
this regime except that it is not obvious that the pic- 
ture of an "adiabatic" continuation from the gapped in- 
sulating state immediately to a gapped superconducting 
state would be valid in a real sample. We will attempt 
to address this issue from another direction by studying 
a heterostructure geometry in Section IIVCI in which we 
can generate a convergent, vortex lattice solution by in 
homogeneously doping the system i.e. high-doping on 
the surface and low-doping in the bulk. 

In the second distinct regime, once is tuned beyond 
|/io|, then significant s-wave pairing begins to develop. 
Because of the particle-hole constraint of the BdG quasi- 
particle spectrum, the energies appear in ±E pairs. The 
lowest energy branches are nearly 2 x x Ny fold degen- 
erate. This degeneracy clearly indicates that these states 
are in-gap vortex states as there is essentially one for each 
vortex. As the chemical potential is pushed more into the 
valence band, the lowest energy branch approaches zero 
energy and at critical chemical potential the par- 
ticle and hole branches cross indicating the location of 
the vortex phase transition. In the weak-pairing treat- 
ment the critical chemical potential is independent of 
the value of the attractive potential \U\ and if we re- 
peat their analysis for our choice of parameters, we find 
a weak-pairing estimate of = 1.35 eV. In our case, 
as the interaction strength is quite large we are not in 
the weak-pairing limit and the critical chemical poten- 
tial depends on the attractive potential. At \U\ =2 eV, 
2.5 eV and 2.8 eV, 1.26 eV, 1.22 eV and 1.2 eV, 

respectively. A stronger \U\ gives a smaller value of 
and it approaches to the weak paring limit as we decrease 
the magnitude of interaction. Since the phenomenon sur- 
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vives the weak-pairing limit it is possible that the vortex 
topological phase transition could also be observed in a 
strong-pairing atomic limit which is realizable in ultra- 
cold optical lattices^^. 

In Fig. m we show the self-consistent vortex profiles in 
a single unit cell as a function of for \U\ = 2.8 eV. It 
is evident that in all cases around the vortex cores, the 
pairing order parameters are suppressed. Away from the 
vortex cores the pairing order parameters are restored 
to A A = 0.21, 0.213 and 0.29 which are roughly equal 
to the corresponding bulk values at = 0.96 eV, 0.98 
eV and 1.3 eV respectively (c.f. Fig. [T]). In the bulk 
superconducting TI, at larger stronger Cooper pair- 
ing is induced, and the strong superconductivity leads to 
a shorter coherence length (^o = ■^-),^^ and thus a 
smaller vortex size. Therefore, in Fig. lU^b) and (c), we 
see flatter order parameter profiles. However, we find un- 
usual behavior in these figures associated with chemical 
potentials of =0.98 eV and 1.3 eV. At these chemical 
potentials, the vortices do not seem to be as well formed 
as they are when = 0.96 eV. This is due to the nu- 
merical discreteness in our simulations. In our system, 
because of the non-trivial order parameter winding due 
to the vortex, there must be a place where the order pa- 
rameter magnitude vanishes. One can see that in Fig. 
Hlthis only happens for = 0.96. What is happening 
is that the vortex core moves from being centered at a 
lattice vertex to the interior of a plaquette. The order 
parameter then vanishes in the plaquette interior (which 
of course is not seen on our discrete lattice spatial sam- 
pling). In fact, we find that at particular values of the 
chemical potential it becomes energetically favorable for 
the vortex to move its core off of a lattice vertex and into 
the center of a plaquette. This is seen in the energy spec- 
tra in Fig. [3] where we a kink in the spectrum a appears 
where this vortex shift occurs, namely around = 0.97 
eV. We believe this is simply an artifact of our numerical 
technique and does not represent any real physics. 



B. Open Vortex Lines in Vortex Lattices 

Next we turn to the case with open boundary condi- 
tions along the z-direction such that the vortices termi- 
nate on the surfaces. This setting is directly relevant for 
possible experiments where the Majorana vortex modes 
are present at the end of vortex lines. Due to the open 
boundaries, the self-consistent BdG calculations must be 
performed in three spatial dimensions as we cannot ex- 
ploit any translation symmetry in z direction. The mag- 
netic unit cells are x ly = S x 16 sized with /^-layers, 
(usually Iz = 6) and there are Nx x Ny = 40 x 20 mag- 
netic unit cells chosen so that we are simulating 1600 
vortices in the vortex lattice. In this section, we only 
consider \U\ =2.8 eV and > |/io| where the supercon- 
ducting gap is formed and the value of \fio\ is estimated 
from the periodic boundary condition case. Here, we self- 
consistently determine the BdG quasi-particle spectrum 




(c) H=1.3 eV 



0.30 - 




0.27- 




FIG. 4. (color online). The spatial pairing order parameter 
profiles |AA,f^| within a unit cell at \U\ = 2.8 eV and at (a) 
= 0.96 eV (b) \fi\ = 0.98 eV and (c) = 1.3 eV. The ver- 
tical axis represents the pairing order parameter magnitudes 
and the horizontal plane is xy plane. Note the variation of 
the pairing order parameter magnitudes at the unit cell center 

Acenter = AA,reunit cell center for different fll (a) Acenter — 0, 
(b) Acenter — 0-12, and (c) Acenter — 0-28. Acenter — indi- 
cates that the vortex core stays on-site, whereas Acenter 7^ 
means that the vortex core moves off the lattice vertices and 
into the plaquette. We note that the profiles for lA^^f^-] look 
similar. 



in the vortex lattice state For the square vortex 

lattice with Nx x Ny magnetic unit cells, there are NxNy 
magnetic Bloch wavevectors k analogous to the wavevec- 
tors in the Brillouin zone of a xA^^^ x Ny square lattice. 

In Fig. O we present the dispersion of vortex modes 
at four different chemical potentials as has been done 
previously for s-wave and d-wave superconductors^*^. 
The high symmetry points of the square lattice are at 
r = (0,0), X = (tt, 0) and M = (7r,7r), as indicated in 
the inset of Fig. [5jc). There are four low-energy "Majo- 
rana" modes at each momentum which are contributed 
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FIG. 5. (color online). The quasiparticle band structure for 
open vortex lines in the bulk superconducting TI. The inter- 
action strength is chosen at \U\ = 2.8 eV and the chemical 
potentials are (a) = 0.6 eV, (b) = 0.9 eV, (c) \fi\ = 1.0 
eV and (d) = 1.3 eV. The inset in (c) denotes the magnetic 
Brillouin zone for the square vortex lattice. The magnetic unit 
cell sizes are Ix x ly y< Iz = Sa x 16a x 6a. 



by the two- vortices per cell and the two ends of each vor- 
tex hne. For a single magnetic unit cell we would thus 
expect to see one Majorana mode on the two ends of 
each of the two vortex lines giving rise to a total of four 
vortices per cell. In this context, we put the word Ma- 
jorana in quotes because, strictly speaking, the low-lying 
energy states only have true Majorana character if they 
are strictly at zero-energy. In Fig. [5fa), we study the 
quasiparticle bands at =0.6 eV and we find that the 
vortex modes are clearly dispersing. Although the super- 
conducting gap is formed, it remains small and the vortex 
modes of different vortices on the same surface can tun- 
nel laterally and hybridize which leads to the dispersion 
of the vortex core states. As we increase the doping level, 
the lowest energy quasiparticle band flattens as is clear 
in the dispersion plot for = 0.9 eV in Fig. [5fb). This 
happens because of the increasing bulk superconducting 
gap, indicated in Fig. ^c). The vortex core size shrinks 



FIG. 6. (color online). The energy splitting 6E vs thickness 
Iz. The magnetic unit cells are 8 x 16 x at = 1 eV and 
1^1 = 1.1 eV. The Hubbard interaction is \U\ = 2.8 eV. 



which leads to smaller overlap of the modes localized in 
different vortices and suppresses the quasiparticle disper- 
sion. This effect {i.e. increase of superconducting gap by 
increasing the doping) stabilizes the Majorana modes. 
The low-energy, flat quasi-particle bands contain AN^Ny 
nearly-degenerate states coming from the 2NxNy vortex 
Majorana modes on the two distinct surfaces. 

In Fig. [5l we see two clear gap-like behaviors. One 
type in Fig. [5ja) shows gaps at low-energy but with 
strong dispersion while Fig. [5fc) and (d) show clear gaps 
but with flat dispersion. For the fiat-dispersing cases we 
studied the dependence of the energy splitting, SE, on 
the sample thickness. An exponential dependence would 
indicate that the dispersionless gap is a result of the hy- 
bridization of the modes at the end of the vortices be- 
tween two surfaces. Fig. [6] shows the energy splitting SE 
has an exponential decreasing relation with the thickness 
L described as^- 



6E (xe ^- 



(15) 



where denotes the characteristic decay length for the 
Majorana modes. A smaller means more localized 
Majorana bound states. By linear fitting from Fig. [6l 
the characteristic length at = 1 eV is ^rn — 5.46a 
and at = 1.1 eV is — 9.49a. This suggests that 
the Majorana modes are still exponentially localized on 
the surface even though there is a gap in our finite-size 
numerics. Therefore, although the Majorana modes may 
tunnel to the opposite surface, in the thermodynamic 
limit {Iz oo), the Majorana modes are still be bound 
to the surface. Although we do not show it here, we note 
that this is not the case for = 1.3 eV where the gap 
does not decrease exponentially, which is expected since 
this is the trivial regime where it should have a power-law 
decay with inverse thickness due to finite-size splitting. 

The nature of the Majorana modes may be further 
illustrated by studying real space probability distribu- 
tions of the in-gap modes. Fig. [7fa)-(d) depict side-view 
spatial slices of the probability density for the lowest en- 
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FIG. 7. (color online). Spatial slices (side-views in the yz plane at x = =b2a) of the probability density for the Majorana 
modes and order parameter density in the bulk superconducting TI at different /i. Upper panels: (a)-(d) shows the evolution 
of the Majorana mode distributions. Brighter regions represent the higher probability density. Lower panels: (e)-(h) shows the 
distribution of pairing order parameters (|Aa + Ab|). The chemical potentials are = 0.6 eV in (a) and (e), = 0.9 eV in 
(b) and (f), = 1 eV in (c) and (g), and = 1.3 eV in (d) and (h), respectively. 



ergy modes and Figs. [71(e)- (h) show the order param- 
eter distributions in real space. The plots are cut on 
the yz surface at x = ±2a, where the vortex cores are 
approximately located. The Majorana modes (indicated 
by bright regions) are observed and localized around the 
vortex cores close to the surfaces in Figs. [TJa) and (b). 
However, the Majorana mode in Fig. [71(a) spreads more 
widely along the surface than that in Fig. [Tfb). This 
shows that at = 0.6 eV neighboring vortices have 
larger overlap than that at =0.9 eV, which corrobo- 
rates with our quasiparticle spectra that indicate stronger 
dispersion for the former case as shown in Fig. [51(a) due 
to the intra-surface hybridization resulting from the in- 
creased lateral overlap of the Majorana modes. It is also 
interesting to see that around =0.9 eV, the mini-gap 
size of the vortex lines is maximum, [see Fig. [31(c)] which 
is where and Fig. [71(f) shows strong, straight-line vortex 
structures. 

At first, further increases in the chemical potential flat- 
tens the dispersion and strengthens the localization of 
the Majorana modes. However, further increases in the 
chemical potential lead to another tunneling mechanism 
for the Majorana modes. As we have already shown for 
periodic vortex rings (e.g see Fig. [3Kc) as > 1) the 
mini-gap of the vortex core states along the vortex line 
eventually begins to decrease as the critical point is ap- 
proached. For open-boundary conditions this leads to 
increased inter-surface hybridization of the modes at the 
two ends of the vortex lines. This results in the formation 
of gaps due to Majorana mode annihilation on opposite 
surfaces (for thin samples) and in Fig. [51(c) we can see 
that there exists a 5E splitting in the Majorana modes. 



As mentioned and shown in Fig. [6l 3E decreases expo- 
nentially in the thickness of the sample. An important 
feature to note is that, as is clear from the dispersion for 
= 0.9 eV, even though the gap increases the band- 
width of quasiparticles decreases and gets flatter. This 
is an indication that intra-surface tunneling is weaken- 
ing (no 2D hopping on the same surface) and that inter- 
surface tunneling is becoming stronger. We can see this 
in Fig. [TKc), where the Majorana bound states begin to 
leak to the opposite surface. Furthermore, in Fig. [71(d), 
in which case the system is topologically trivial, the low- 
est energy modes, which are no longer Majorana in na- 
ture and gapped by the vortex mini-gap of order A^//i, 
completely penetrate through the bulk at = 1.3 eV 
and lie along the vortex lines. 



Superconductor- Topological Insulator 
Heterostructure 



As mentioned in the introduction, another method to 
realize the Majorana modes is through the proximity ef- 
fect of a topological insulator and an s-wave supercon- 
ductor. Using our self-consistent BdG method, we can 
also study this geometry. By modeling such a structure 
by an inhomogenous doping level we can directly address 
the effect of the penetration of superconducting gap into 
the bulk and the self-consistent formation of Majorana 
bound states in vortices. We imagine a similar proximity- 
induced superconductivity as was first suggested by Fu 
and Kane^^. We choose an inhomogenous system where 
/i(r) in Eq. ([TT]) is layer-dependent. We choose the sur- 
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FIG. 8. (color online). The spatial side view of the pairing 
order parameter distribution of a six- layer the s-wave/TI het- 
erostructure. The interaction strength is chosen at \U\ = 2.8 
eV, and the surface and bulk chemical potentials are = 1 
eV and = 0.55 eV, respectively. The dark blue tubes in- 
dicate the region that pairing is suppressed and form vortex 
lines. The magnetic unit cell sizes are 8a x 16a x 6a 

face chemical potential /i(r) = /i^ = — 1 eV, and the 
chemical potential /i(r) = fiB = —0.55 eV in other bulk 
layers so that both the bulk and surface have non-zero 
density of states. We investigate the case where super- 
conductivity is induced primarily on the surfaces by turn- 
ing on the same attractive interaction across the entire 
sample. The inhomogeneous /i(r) will generate a much 
stronger order parameter on the surfaces than in the bulk 
due to the large difference in chemical potentials. Again 
we choose a uniform magnetic field along the z-direction 
which generates the vortex lattice. 

There is another reason to consider this system be- 
yond simply the presence of Major ana fermions. One of 
the major obstacles in our bulk calculation is the non- 
convergence of a stable vortex solution when the order 
parameter magnitude is very small. As mentioned, when 
l/i^l < < \fio\ despite the presence of gapless elec- 
trons at the Fermi-level, the density of states is not large 
enough to form a sizable superconducting gap (at least 
for system sizes we consider) and the doped TI remains a 
gapless metal. We can counter-act this problem by using 
the superconductor-TI heterostructure geometry which 
acts to pin the vortices with strong superconductivity at 
the surface (highly-doped region) thereby stabilizing the 
solution. Fig. [8] shows the spatial side view of the result- 
ing self-consistent order parameter profile of a six-layer 
heterostructure. We see no evidence of superconductiv- 
ity in the bulk and roughly uniform superconductivity in 
the surface which is interrupted near the vortices. The 
resulting calculation in Fig. [9)3 shows that Majorana sur- 
face modes still remain even though the superconduct- 
ing order parameter in the bulk is exponentially small 
compared to the surface. The six-layer heterostructure 
can roughly approximate the case of two vortices exist- 
ing in a four-layer bulk TI which is uniformly doped with 
= 0.55 eV. This was a region of interest that we could 
not access in our bulk calculation due to finite-size com- 
plications and which we can, admittedly only roughly. 
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FIG. 9. (color online), (a) The quasiparticle band spectrum 
for six- layer the s-wave/TI heterostructure. (b) The spatial 
side view of the Majorana modes. The interaction strength is 
chosen at \U\ = 2.8 eV and = 1 eV and l/i^l = 0.55 eV. 



learn about by stabilizing the vortex solution using higher 
surface doping. 

In Fig. ini^a), we present the quasi-particle band spec- 
trum for the superconductor-TI heterostructure with a 
square vortex lattice. Within the superconducting gap, 
there exist two prominent low-energy modes which are 
doubly degenerate whose energies are split away from 
zero energy. Although not shown here, we find that the 
energy splitting 6E also has an exponential decay with 
increasing sample thickness Iz. This indicates that the 
low-energy modes are exponentially localized at the su- 
perconducting surface and the Majorana fermions can 
stably reside at the surfaces in the thermodynamic limit, 
i.e. Iz ^ oo. This is indicative that that the vortices in 
the low-doping regime (|/i^| < < |//o|) also support 
Majorana fermions in the bulk superconducting TI. In 
Fig. [njb) we show the probability density of the low en- 
ergy modes and see the tight localization of the resulting 
modes on the surface as one would expect for Majorana 
modes formed within a well- formed superconducting gap. 
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FIG. 10. (color online) The comparison between the midgap 

sizes Um 

(T) and ksT. The intersection occurs at T = To ^ 
0.025 K indicating that as T < To, the Majorana modes can 
stably exist on the surface of the doped topological insulators. 



exist on the surfaces and can be detected experimentally. 
The numerical result is shown in Fig. [TOl From Ref. |52|, 
in Bi2Se3, 6eF ^ 0.25 eV which results in an estimate for 
the critical temperature for the observation of Majorana 
modes to be Tq ~ 0.025 K. Therefore, we can provide 
a rough estimate that at T < 0.025 K, the Majorana 
modes can stably exist on the surface of the doped topo- 
logical insulators and may be detectable experimentally. 
This number is quite small and indicates one would need 
to optimize materials properties in order to hope for ob- 
servation. The results for Heusler materials or materials 
with similar electronic structure to bulk HgTe may pro- 
vide more promising alternatives^^ due to the differences 
in the sustainable levels of doping. 



VI. CONCLUSION 



V. VORTEX MAJORANA MODES AT FINITE 
TEMPERATURES 

In our previous analysis, we have neglected the role of 
temperature in our analysis. In this section, we provide 
a rough estimate of the temperature at which one could 
observe the Majorana fermions experimentally. In the 
BCS theory the temperature dependence of the gap size 
is determined through)^ 



1 



N{0)V 



hujc 



tanh 



2kBT 



d^. (16) 



Here A^(0) is the density of states at the Fermi- level, 
U is the interaction coupling, and Uc is the Debye fre- 
quency. The finite temperature gap A(T) may only be 
determined numerically. The combination of N{0)U can 
be estimated from the critical temperature, T^, and the 
Debye frequency, cjc, for Cu-doped Bi2Se3 via: 



N{0)U 



In 



(17) 



For Bi2Se3, the critical temperature is Tc = 3.8K,^^^^'^ 
and the Debye temperature huOc/kB is 180 K^^. With 
N{0)U determined from Eq. (pT|) . one can calculate the 
temperature dependence of the gap numerically using Eq. 
(p!6|) . To observe the Majorana modes at finite tempera- 
ture, the mini-gap size of the vortex lines Sm{T) should be 
stable against the thermal fluctuations: Sm{T) < ksT. 
The temperature at which this occurs may be estimated 
as 



To 



5m{T) 7rA{Ty 



(18) 



where Sep = |/i — /i5| where /i^ denotes the bulk band 
edge. When T < Tq^ the Majorana modes can stably 



In summary, we performed self-consistent Bogoliubov- 
de-Gennes calculations to study properties of vortices in 
doped topological insulators that become superconduct- 
ing. Through the use of our numerics, we studied the 
physics of Majorana fermions in vortex lattices beyond 
the strict weak-coupling hmit, and the resulting vortex 
phase transitions between a topological and trivial state. 
We have shown that the quasi-particle band spectra offers 
evidence that there exists an optimal regime in chemical 
potential where the Majorana fermions can stably reside 
even in a finite thickness system. There also exists other 
regimes where the Majorana fermions do not stably ex- 
ist on the system surfaces because of intra- and inter- 
surface hybridization between the vortex modes. Fur- 
thermore, we also showed that, through the use of the 
analogous s-wave-TI heterostructure, that TIs with bulk 
superconductivity containing finite carrier density but in- 
sufficient superconducting pairing strength can host Ma- 
jorana fermions on the surface. Similar to the bulk super- 
conducting case, the Majorana modes can also leak into 
the bulk and annihilate with the other surface. However, 
the tunneling behavior exhibits the usual exponential de- 
cay with thickness and we conclude that the Majorana 
fermions can survive for thick samples. Unfortunately, 
the simple estimates we made for a viable temperature 
range in which Majorana modes may be observed indi- 
cate that superconducting Cu-Bi2Se3 may not provide a 
good candidate even if the doping level can be tuned to 
the topological vortex phase. 
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